Based on: Handbook of Parametric and Nonparametric Statistical Procedures, 5th edition, Chapter 11

using Random
Random.seed!(783376894512)
using Distributions
using Intervals
using Plots
using LinearAlgebra
n_sim = 1_000_000 # number of simulations to verify theoretical distributions

True model

n1 = 4
n2 = 5
σₜ = 2.0
μₜ1 = 3.0
μₜ2 = 4.0
Mₜ1 = Product([Normal(μₜ1, σₜ) for _ in 1:n1]) # assumed known to be Normal and i.i.d.
Mₜ2 = Product([Normal(μₜ2, σₜ) for _ in 1:n2]) # assumed known to be Normal and i.i.d.
generate_data(M1, M2) = (rand(M1), rand(M2))

Hypothesis testing

function test_statistic(y1, y2, Δμ₀)
    n1 = length(y1)
    n2 = length(y2)
    μ̂1 = sum(y1) / n1
    μ̂2 = sum(y2) / n2
    σ̂²1 = sum((y1 .- μ̂1) .^ 2) / (n1 - 1)
    σ̂²2 = sum((y2 .- μ̂2) .^ 2) / (n2 - 1)
    σ̂²p = ((n1 - 1) * σ̂²1 + (n2 - 1) * σ̂²2) / (n1 + n2 - 2)
    return (μ̂1 - μ̂2 - Δμ₀) / (sqrt(σ̂²p) * sqrt(1 / n1 + 1 / n2))
end

Null Hypothesis is true

Δμ₀ = -1.0
Δμ₀ == μₜ1 - μₜ2

all(mean(Mₜ1) .== μₜ1)
all(mean(Mₜ2) .== μₜ2)
all(sqrt.(var((Mₜ1))) .== σₜ)
all(sqrt.(var((Mₜ2))) .== σₜ)

simulated_statistics = [test_statistic(generate_data(Mₜ1, Mₜ2)..., Δμ₀) for _ in 1:n_sim]
distribution_under_H₀ = t -> pdf(TDist(n1 + n2 - 2), t)

Code

histogram(simulated_statistics; normalize=true, label="Simulated")
tmin = mean(simulated_statistics) - 3 * std(simulated_statistics)
tmax = mean(simulated_statistics) + 3 * std(simulated_statistics)
t_plot = tmin:((tmax - tmin) / 1000):tmax
plot!(distribution_under_H₀, t_plot; lw=5, label="Theoretical")
plot!(;
    title="Distribution of two independent samples T test statistic, \n when Null Hypothesis is correct",
    xlimits=(tmin, tmax),
    xlabel="t",
    ylabel="pdf",
    grid=false,
    legendfontsize=12,
    tickfontsize=12,
    guidefontsize=16,
)

Null Hypothesis is true, but assumptions don’t hold

Observations are non-normal

M1_non_normal = Product([
    MixtureModel(
        [Normal(μₜ1 + sqrt(3.99), 0.1), Normal(μₜ1 - sqrt(3.99), 0.1)], [0.5, 0.5]
    ) for _ in 1:n1
])
M2_non_normal = Product([
    MixtureModel(
        [Normal(μₜ2 + sqrt(3.99), 0.1), Normal(μₜ2 - sqrt(3.99), 0.1)], [0.5, 0.5]
    ) for _ in 1:n2
])
all(mean(M1_non_normal) .== μₜ1)
all(mean(M2_non_normal) .== μₜ2)
all(sqrt.(var(M1_non_normal)) .== σₜ)
all(sqrt.(var(M2_non_normal)) .== σₜ)
simulated_statistics = [
    test_statistic(generate_data(M1_non_normal, M2_non_normal)..., Δμ₀) for _ in 1:n_sim
]

Code

histogram(simulated_statistics; normalize=true, label="Simulated")
tmin = mean(simulated_statistics) - 3 * std(simulated_statistics)
tmax = mean(simulated_statistics) + 3 * std(simulated_statistics)
t_plot = tmin:((tmax - tmin) / 1000):tmax
plot!(distribution_under_H₀, t_plot; lw=5, label="Theoretical")
plot!(;
    title="Distribution of two independent samples T test statistic, \n under H₀, but non-normal observations",
    xlimits=(tmin, tmax),
    xlabel="t",
    ylabel="pdf",
    grid=false,
    legendfontsize=12,
    tickfontsize=12,
    guidefontsize=16,
)

Observations are non-independent

M1_non_independent = MvNormal(
    [μₜ1, μₜ1, μₜ1, μₜ1, μₜ1],
    [
        σₜ^2 0.5*σₜ^2 0 0 0
        0.5*σₜ^2 σₜ^2 0 0 0
        0 0 σₜ^2 0 0
        0 0 0 σₜ^2 0
        0 0 0 0 σₜ^2
    ],
)
Σ2 = collect(Diagonal(fill(σₜ^2, n2)))
Σ2[2, 1] = Σ2[1, 2] = 0.5 * σₜ^2
M2_non_independent = MvNormal(fill(μₜ2, n2), Σ2)

all(mean(M1_non_independent) .== μₜ1)
all(mean(M2_non_independent) .== μₜ2)
all(sqrt.(var(M1_non_independent)) .== σₜ)
all(sqrt.(var(M2_non_independent)) .== σₜ)
simulated_statistics = [
    test_statistic(generate_data(M2_non_independent, M2_non_independent)..., Δμ₀) for
    _ in 1:n_sim
]

Code

histogram(simulated_statistics; normalize=true, label="Simulated")
tmin = mean(simulated_statistics) - 3 * std(simulated_statistics)
tmax = mean(simulated_statistics) + 3 * std(simulated_statistics)
t_plot = tmin:((tmax - tmin) / 1000):tmax
plot!(distribution_under_H₀, t_plot; lw=5, label="Theoretical")
plot!(;
    title="Distribution of two independent samples T test statistic, \n under H₀, but non-independent observations",
    xlimits=(tmin, tmax),
    xlabel="t",
    ylabel="pdf",
    grid=false,
    legendfontsize=12,
    tickfontsize=12,
    guidefontsize=16,
)

Observations are non-identical

M1_non_identical = Product([Normal(μₜ1, i * σₜ) for i in 1:n1])
M2_non_identical = Product([Normal(μₜ2, i * σₜ) for i in 1:n1]) # note the n1 here to differentiate from the next counterexample
all(mean(M1_non_identical) .== μₜ1)
all(mean(M2_non_identical) .== μₜ2)
any(sqrt.(var(M1_non_identical)) .== σₜ)
any(sqrt.(var(M2_non_identical)) .== σₜ)

simulated_statistics = [
    test_statistic(generate_data(M1_non_identical, M2_non_identical)..., Δμ₀) for
    _ in 1:n_sim
]

Code

histogram(simulated_statistics; normalize=true, label="Simulated")
tmin = mean(simulated_statistics) - 3 * std(simulated_statistics)
tmax = mean(simulated_statistics) + 3 * std(simulated_statistics)
t_plot = tmin:((tmax - tmin) / 1000):tmax
plot!(distribution_under_H₀, t_plot; lw=5, label="Theoretical")
plot!(;
    title="Distribution of two independent samples T test statistic, \n under H₀, but non-identical observations",
    xlimits=(tmin, tmax),
    xlabel="t",
    ylabel="pdf",
    grid=false,
    legendfontsize=12,
    tickfontsize=12,
    guidefontsize=16,
)

Two samples have a different variance

M2_different_variance = Product([Normal(μₜ1, 2 * σₜ) for _ in 1:n2])
all(mean(M2_different_variance) .== μₜ1)
all(sqrt.(var(M2_different_variance)) .!= σₜ)

simulated_statistics = [
    test_statistic(generate_data(Mₜ1, M2_different_variance)..., Δμ₀) for _ in 1:n_sim
]

Code

histogram(simulated_statistics; normalize=true, label="Simulated")
tmin = mean(simulated_statistics) - 3 * std(simulated_statistics)
tmax = mean(simulated_statistics) + 3 * std(simulated_statistics)
t_plot = tmin:((tmax - tmin) / 1000):tmax
plot!(distribution_under_H₀, t_plot; lw=5, label="Theoretical")
plot!(;
    title="Distribution of two independent samples T test statistic, \n under H₀, but different sample variances",
    xlimits=(tmin, tmax),
    xlabel="t",
    ylabel="pdf",
    grid=false,
    legendfontsize=12,
    tickfontsize=12,
    guidefontsize=16,
)

Null Hypothesis is false

Δμ₀ = 2.0
Δμ₀ != μₜ1 - μₜ2

simulated_statistics = [test_statistic(generate_data(Mₜ1, Mₜ2)..., Δμ₀) for _ in 1:n_sim]
distribution_under_H₁ =
    tnc -> pdf(NoncentralT(n1 + n2 - 2, (μₜ1 - μₜ2 - Δμ₀) / √(σₜ^2 / n1 + σₜ^2 / n2)), tnc)

Code

histogram(simulated_statistics; normalize=true, label="Simulated")
tmin = mean(simulated_statistics) - 3 * std(simulated_statistics)
tmax = mean(simulated_statistics) + 3 * std(simulated_statistics)
t_plot = tmin:((tmax - tmin) / 1000):tmax
plot!(distribution_under_H₁, t_plot; lw=5, label="Theoretical")
plot!(;
    title="Distribution of two independent samples T test statistic, \n when Null Hypothesis is incorrect",
    xlimits=(tmin, tmax),
    xlabel="t",
    ylabel="pdf",
    grid=false,
    legendfontsize=12,
    tickfontsize=12,
    guidefontsize=16,
)

Null Hypothesis is false, but assumptions don’t hold

Same counter examples as above.

Code

simulated_statistics = [
    test_statistic(generate_data(M1_non_normal, M2_non_normal)..., Δμ₀) for _ in 1:n_sim
]

1000000-element Vector{Float64}:
 -2.7319646451633637
 -1.1635565113397672
 -2.5973138575913843
 -3.1753393673498325
 -0.6945610789746434
 -0.5825206489634959
 -0.9541786713792212
 -0.6833793462505765
 -1.6870049418296202
 -1.1047206335419304
  ⋮
 -1.0827951572184749
 -2.4820903467430555
 -1.342302654041741
 -3.0872701991748603
 -1.0679875565707986
 -3.2196581992897255
 -0.7061209569844831
 -1.6155670827847624
 -0.580786774694805

Code

histogram(simulated_statistics; normalize=true, label="Simulated")
tmin = mean(simulated_statistics) - 3 * std(simulated_statistics)
tmax = mean(simulated_statistics) + 3 * std(simulated_statistics)
t_plot = tmin:((tmax - tmin) / 1000):tmax
plot!(distribution_under_H₁, t_plot; lw=5, label="Theoretical")
plot!(;
    title="Distribution of two independent samples T test statistic, \n under H₁, but non-normal observations",
    xlimits=(tmin, tmax),
    xlabel="t",
    ylabel="pdf",
    grid=false,
    legendfontsize=12,
    tickfontsize=12,
    guidefontsize=16,
)

Code

simulated_statistics = [
    test_statistic(generate_data(M2_non_independent, M2_non_independent)..., Δμ₀) for
    _ in 1:n_sim
]

histogram(simulated_statistics; normalize=true, label="Simulated")
tmin = mean(simulated_statistics) - 3 * std(simulated_statistics)
tmax = mean(simulated_statistics) + 3 * std(simulated_statistics)
t_plot = tmin:((tmax - tmin) / 1000):tmax
plot!(distribution_under_H₁, t_plot; lw=5, label="Theoretical")
plot!(;
    title="Distribution of two independent samples T test statistic, \n under H₁, but non-independent observations",
    xlimits=(tmin, tmax),
    xlabel="t",
    ylabel="pdf",
    grid=false,
    legendfontsize=12,
    tickfontsize=12,
    guidefontsize=16,
)

Code

simulated_statistics = [
    test_statistic(generate_data(M1_non_identical, M2_non_identical)..., Δμ₀) for
    _ in 1:n_sim
]

histogram(simulated_statistics; normalize=true, label="Simulated")
tmin = mean(simulated_statistics) - 3 * std(simulated_statistics)
tmax = mean(simulated_statistics) + 3 * std(simulated_statistics)
t_plot = tmin:((tmax - tmin) / 1000):tmax
plot!(distribution_under_H₁, t_plot; lw=5, label="Theoretical")
plot!(;
    title="Distribution of two independent samples T test statistic, \n under H₁, but non-identical observations",
    xlimits=(tmin, tmax),
    xlabel="t",
    ylabel="pdf",
    grid=false,
    legendfontsize=12,
    tickfontsize=12,
    guidefontsize=16,
)

Code

simulated_statistics = [
    test_statistic(generate_data(Mₜ1, M2_different_variance)..., Δμ₀) for _ in 1:n_sim
]

histogram(simulated_statistics; normalize=true, label="Simulated")
tmin = mean(simulated_statistics) - 3 * std(simulated_statistics)
tmax = mean(simulated_statistics) + 3 * std(simulated_statistics)
t_plot = tmin:((tmax - tmin) / 1000):tmax
plot!(distribution_under_H₁, t_plot; lw=5, label="Theoretical")
plot!(;
    title="Distribution of two independent samples T test statistic, \n under H₁, but different sample variances",
    xlimits=(tmin, tmax),
    xlabel="t",
    ylabel="pdf",
    grid=false,
    legendfontsize=12,
    tickfontsize=12,
    guidefontsize=16,
)

Confidence Interval

α = 0.05 # proportion of CIs which should fail to cover true values
function confidence_interval(y1, y2, α)
    n1 = length(y1)
    n2 = length(y2)
    μ̂1 = sum(y1) / n1
    μ̂2 = sum(y2) / n2
    σ̂²1 = sum((y1 .- μ̂1) .^ 2) / (n1 - 1)
    σ̂²2 = sum((y2 .- μ̂2) .^ 2) / (n2 - 1)
    σ̂²p = ((n1 - 1) * σ̂²1 + (n2 - 1) * σ̂²2) / (n1 + n2 - 2)
    L = μ̂1 - μ̂2 + quantile(TDist(n1 + n2 - 2), α / 2) * (√(σ̂²p * (1 / n1 + 1 / n2)))
    U = μ̂1 - μ̂2 + quantile(TDist(n1 + n2 - 2), 1 - α / 2) * (√(σ̂²p * (1 / n1 + 1 / n2)))
    return Interval(L, U)
end

simulated_CIs = [confidence_interval(generate_data(Mₜ1, Mₜ2)..., α) for _ in 1:n_sim]
coverage = count([μₜ1 - μₜ2 in CI for CI in simulated_CIs]) / n_sim

Code

bar(["inside CI", "outside CI"], [coverage, 1 - coverage])
plot!(;
    title="Simulated Coverage of 95% confidence interval \n corresponding to two independent samples T test",
    ylimits=(0.0, 1.0),
    yticks=[0.05, 0.95],
    ylabel="proportion",
    grid=true,
    legend=false,
    legendfontsize=12,
    tickfontsize=12,
    guidefontsize=16,
)

assumptions don’t hold

Same counter examples as above.

Code

simulated_CIs = [confidence_interval(generate_data(M1_non_normal, M2_non_normal)..., α) for _ in 1:n_sim]
coverage = count([μₜ1 - μₜ2 in CI for CI in simulated_CIs]) / n_sim
bar(["inside CI", "outside CI"], [coverage, 1 - coverage])
plot!(;
    title="Simulated Coverage of 95% confidence interval \n corresponding to two independent samples T test \n but non-normal observations",
    ylimits=(0.0, 1.0),
    yticks=[0.05, 0.95],
    ylabel="proportion",
    grid=true,
    legend=false,
    legendfontsize=12,
    tickfontsize=12,
    guidefontsize=16,
    #top_margin = 25Plots.px,
)

Code

simulated_CIs = [confidence_interval(generate_data(M1_non_independent, M2_non_independent)..., α) for _ in 1:n_sim]
coverage = count([μₜ1 - μₜ2 in CI for CI in simulated_CIs]) / n_sim

bar(["inside CI", "outside CI"], [coverage, 1 - coverage])
plot!(;
    title="Simulated Coverage of 95% confidence interval \n corresponding to two independent samples T test \n but non-independent observations",
    ylimits=(0.0, 1.0),
    yticks=[0.05, 0.95],
    ylabel="proportion",
    grid=true,
    legend=false,
    legendfontsize=12,
    tickfontsize=12,
    guidefontsize=16,
)

Code

simulated_CIs = [confidence_interval(generate_data(M1_non_identical, M2_non_identical)..., α) for _ in 1:n_sim]
coverage = count([μₜ1 - μₜ2 in CI for CI in simulated_CIs]) / n_sim

bar(["inside CI", "outside CI"], [coverage, 1 - coverage])
plot!(;
    title="Simulated Coverage of 95% confidence interval \n corresponding to two independent samples T test \n but non-identical observations",
    ylimits=(0.0, 1.0),
    yticks=[0.05, 0.95],
    ylabel="proportion",
    grid=true,
    legend=false,
    legendfontsize=12,
    tickfontsize=12,
    guidefontsize=16,
)

Code

simulated_CIs = [confidence_interval(generate_data(Mₜ1, M2_different_variance)..., α) for _ in 1:n_sim]
coverage = count([μₜ1 - μₜ2 in CI for CI in simulated_CIs]) / n_sim

bar(["inside CI", "outside CI"], [coverage, 1 - coverage])
plot!(;
    title="Simulated Coverage of 95% confidence interval \n corresponding to two independent samples T test \n but different sample variances",
    ylimits=(0.0, 1.0),
    yticks=[0.05, 0.95],
    ylabel="proportion",
    grid=true,
    legend=false,
    legendfontsize=12,
    tickfontsize=12,
    guidefontsize=16,
)